##########################################################################################

library(data.table)
library(optparse)
library(phastCons100way.UCSC.hg38)
library(ggplot2)

##########################################################################################

option_list <- list(
    make_option(c("--somatic_peak_file"), type = "character"),
    make_option(c("--somatic_peak_peak_gene_file"), type = "character"),
    make_option(c("--germ_peak_file"), type = "character"),
    make_option(c("--germ_peak_peak_gene_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    somatic_peak_file <- "~/20231121_singleMuti/results/qc_atac_v3/somatic/testis_combined_peak.qc.tsv"
    somatic_peak_peak_gene_file <- "~/20231121_singleMuti/results/celltype_plot/peak2gene/somatic/peakToGeneHeatmap_LabelClust_k9.tsv"
    germ_peak_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ/testis_combined_peak.qc.tsv"
    germ_peak_peak_gene_file <- "~/20231121_singleMuti/results/celltype_plot/peak2gene/germ/peakToGeneHeatmap_LabelClust_k12.tsv"

    out_path  <- "~/20231121_singleMuti/results/celltype_plot/peak2gene"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

somatic_peak_file <- opt$somatic_peak_file
somatic_peak_peak_gene_file <- opt$somatic_peak_peak_gene_file
germ_peak_file <- opt$germ_peak_file
germ_peak_peak_gene_file <- opt$germ_peak_peak_gene_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

somatic_peak <- fread(somatic_peak_file)
somatic_peak_peak_gene <- fread(somatic_peak_peak_gene_file)
germ_peak <- fread(germ_peak_file)
germ_peak_peak_gene <- fread(germ_peak_peak_gene_file)

phast <- phastCons100way.UCSC.hg38

###########################################################################################
## 注释细胞保守型
computeCons <- function( dat_peak = dat_peak , dat_peak_gene = dat_peak_gene , peak_type = peak_type ){
    dat_peak$peak <- paste0(dat_peak$seqnames , ":" , dat_peak$start , "-" , dat_peak$end)
    ## 判断是否为peak2gene
    dat_peak$peak2gene <- ifelse( dat_peak$peak %in% dat_peak_gene$peak , "TRUE" , "FALSE" )

    ## 注释保守性评分
    allPeaksGR <- makeGRangesFromDataFrame(dat_peak , keep.extra.columns=TRUE)
    names(allPeaksGR) <- allPeaksGR$peak
    allPeaksCons <- gscores(phast, allPeaksGR, summaryFun=mean) # Mean is the default summaryFun
    allPeaksCons$peak_type <- peak_type
    allPeaksCons <- data.frame(allPeaksCons)

    return(allPeaksCons)
}

dat_peak <- somatic_peak
dat_peak_gene <- somatic_peak_peak_gene
peak_type <- "somatic"
dat_somatic <- computeCons( dat_peak = dat_peak , dat_peak_gene = dat_peak_gene , peak_type = peak_type )

dat_peak <- germ_peak
dat_peak_gene <- germ_peak_peak_gene
peak_type <- "germ"
dat_germ <- computeCons( dat_peak = dat_peak , dat_peak_gene = dat_peak_gene , peak_type = peak_type )

dat_combine <- rbind( dat_germ , dat_somatic )
colnames(dat_combine)[which(colnames(dat_combine)=="default")]="phastCons100way"

out_file <- paste0( out_path , "/germ_somtic-peak-gene_phastCons100way.tsv" )
write.table( dat_combine , out_file , row.names = F , sep = "\t" , quote = F )

###########################################################################################
## 画图
all_data <- dat_combine
all_data$peak2gene_type <- ifelse( all_data$peak2gene == "TRUE" , "peaklink" , "other" )
all_data$peak2gene_type <- factor(all_data$peak2gene_type, levels = c("peaklink" , "other"))

p1<-ggplot(data = all_data , aes(x = peak_type , y= phastCons100way , fill = peak2gene_type))+
  geom_boxplot(width=0.6,
               position = position_dodge(0.6))+
  theme_bw()+
  theme(panel.grid = element_blank())+
  scale_fill_manual(values = c("yellowgreen", "violetred1"))

out_file <- paste0( out_path , "/germ_somtic-peak-gene_phastCons100way.pdf" )
ggsave( out_file , p1 )
  